<html>
<head>
<title>eigen.c</title>
</head>
<body>
<tt><b><i>//-----------------------------------------------------------------------------</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute&nbsp;eigenvalues&nbsp;and&nbsp;eigenvectors</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Input:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;stack[tos&nbsp;-&nbsp;1]&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;symmetric&nbsp;matrix</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Output:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;D&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;diagnonal&nbsp;matrix</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Q&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;eigenvector&nbsp;matrix</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;D&nbsp;and&nbsp;Q&nbsp;have&nbsp;the&nbsp;property&nbsp;that</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A&nbsp;==&nbsp;dot(transpose(Q),D,Q)</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;where&nbsp;A&nbsp;is&nbsp;the&nbsp;original&nbsp;matrix.</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;The&nbsp;eigenvalues&nbsp;are&nbsp;on&nbsp;the&nbsp;diagonal&nbsp;of&nbsp;D.</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;The&nbsp;eigenvectors&nbsp;are&nbsp;row&nbsp;vectors&nbsp;in&nbsp;Q.</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;The&nbsp;eigenvalue&nbsp;relation</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A&nbsp;X&nbsp;=&nbsp;lambda&nbsp;X</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;can&nbsp;be&nbsp;checked&nbsp;as&nbsp;follows:</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;lambda&nbsp;=&nbsp;D[1,1]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;X&nbsp;=&nbsp;Q[1]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;dot(A,X)&nbsp;-&nbsp;lambda&nbsp;X</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//-----------------------------------------------------------------------------</i></b></tt><br>
<tt></tt><br>
<tt>#include&nbsp;"defs.h"</tt><br>
<tt></tt><br>
<tt>#define&nbsp;D(i,&nbsp;j)&nbsp;yydd[n&nbsp;*&nbsp;(i)&nbsp;+&nbsp;(j)]</tt><br>
<tt>#define&nbsp;Q(i,&nbsp;j)&nbsp;yyqq[n&nbsp;*&nbsp;(i)&nbsp;+&nbsp;(j)]</tt><br>
<tt></tt><br>
<tt>extern&nbsp;void&nbsp;<a href="tensor.c.html#copy_tensor">copy_tensor</a>(void);</tt><br>
<tt>static&nbsp;void&nbsp;eigen(int);</tt><br>
<tt>static&nbsp;int&nbsp;check_arg(void);</tt><br>
<tt>static&nbsp;int&nbsp;step(void);</tt><br>
<tt>static&nbsp;void&nbsp;step2(int,&nbsp;int);</tt><br>
<tt>static&nbsp;int&nbsp;n;</tt><br>
<tt>static&nbsp;double&nbsp;*yydd,&nbsp;*yyqq;</tt><br>
<tt></tt><br>
<tt>void</tt><br>
<tt><a name="eval_eigen">eval_eigen</a>(void)</tt><br>
<tt>{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(check_arg()&nbsp;==&nbsp;0)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="run.c.html#stop">stop</a>("eigen:&nbsp;argument&nbsp;is&nbsp;not&nbsp;a&nbsp;square&nbsp;matrix");</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;eigen(EIGEN);</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p1&nbsp;=&nbsp;<a href="symbol.c.html#usr_symbol">usr_symbol</a>("D");</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="symbol.c.html#set_binding">set_binding</a>(p1,&nbsp;p2);</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p1&nbsp;=&nbsp;<a href="symbol.c.html#usr_symbol">usr_symbol</a>("Q");</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="symbol.c.html#set_binding">set_binding</a>(p1,&nbsp;p3);</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="stack.c.html#push">push</a>(<a href="defs.h.html#symbol">symbol</a>(NIL));</tt><br>
<tt>}</tt><br>
<tt></tt><br>
<tt>void</tt><br>
<tt><a name="eval_eigenval">eval_eigenval</a>(void)</tt><br>
<tt>{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(check_arg()&nbsp;==&nbsp;0)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="symbol.c.html#push_symbol">push_symbol</a>(EIGENVAL);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="stack.c.html#push">push</a>(p1);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="list.c.html#list">list</a>(2);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;return;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;eigen(EIGENVAL);</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="stack.c.html#push">push</a>(p2);</tt><br>
<tt>}</tt><br>
<tt></tt><br>
<tt>void</tt><br>
<tt><a name="eval_eigenvec">eval_eigenvec</a>(void)</tt><br>
<tt>{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(check_arg()&nbsp;==&nbsp;0)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="symbol.c.html#push_symbol">push_symbol</a>(EIGENVEC);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="stack.c.html#push">push</a>(p1);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="list.c.html#list">list</a>(2);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;return;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;eigen(EIGENVEC);</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="stack.c.html#push">push</a>(p3);</tt><br>
<tt>}</tt><br>
<tt></tt><br>
<tt>static&nbsp;int</tt><br>
<tt>check_arg(void)</tt><br>
<tt>{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;int&nbsp;i,&nbsp;j;</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="stack.c.html#push">push</a>(<a href="defs.h.html#cadr">cadr</a>(p1));</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="eval.c.html#eval">eval</a>();</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="float.c.html#yyfloat">yyfloat</a>();</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="eval.c.html#eval">eval</a>();</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p1&nbsp;=&nbsp;<a href="stack.c.html#pop">pop</a>();</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(!<a href="defs.h.html#istensor">istensor</a>(p1))</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;return&nbsp;0;</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(p1-&gt;u.tensor-&gt;ndim&nbsp;!=&nbsp;2&nbsp;||&nbsp;p1-&gt;u.tensor-&gt;dim[0]&nbsp;!=&nbsp;p1-&gt;u.tensor-&gt;dim[1])</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="run.c.html#stop">stop</a>("eigen:&nbsp;argument&nbsp;is&nbsp;not&nbsp;a&nbsp;square&nbsp;matrix");</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n&nbsp;=&nbsp;p1-&gt;u.tensor-&gt;dim[0];</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(i&nbsp;=&nbsp;0;&nbsp;i&nbsp;&lt;&nbsp;n;&nbsp;i++)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(j&nbsp;=&nbsp;0;&nbsp;j&nbsp;&lt;&nbsp;n;&nbsp;j++)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(!<a href="defs.h.html#isdouble">isdouble</a>(p1-&gt;u.tensor-&gt;elem[n&nbsp;*&nbsp;i&nbsp;+&nbsp;j]))</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="run.c.html#stop">stop</a>("eigen:&nbsp;matrix&nbsp;is&nbsp;not&nbsp;numerical");</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(i&nbsp;=&nbsp;0;&nbsp;i&nbsp;&lt;&nbsp;n&nbsp;-&nbsp;1;&nbsp;i++)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(j&nbsp;=&nbsp;i&nbsp;+&nbsp;1;&nbsp;j&nbsp;&lt;&nbsp;n;&nbsp;j++)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(fabs(p1-&gt;u.tensor-&gt;elem[n&nbsp;*&nbsp;i&nbsp;+&nbsp;j]-&gt;u.d&nbsp;-&nbsp;p1-&gt;u.tensor-&gt;elem[n&nbsp;*&nbsp;j&nbsp;+&nbsp;i]-&gt;u.d)&nbsp;&gt;&nbsp;1e-10)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="run.c.html#stop">stop</a>("eigen:&nbsp;matrix&nbsp;is&nbsp;not&nbsp;symmetrical");</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;return&nbsp;1;</tt><br>
<tt>}</tt><br>
<tt></tt><br>
<tt><b><i>//-----------------------------------------------------------------------------</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Input:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p1&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;matrix</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Output:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p2&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;eigenvalues</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p3&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;eigenvectors</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//-----------------------------------------------------------------------------</i></b></tt><br>
<tt></tt><br>
<tt>static&nbsp;void</tt><br>
<tt>eigen(int&nbsp;op)</tt><br>
<tt>{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;int&nbsp;i,&nbsp;j;</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;malloc&nbsp;working&nbsp;vars</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;yydd&nbsp;=&nbsp;(double&nbsp;*)&nbsp;malloc(n&nbsp;*&nbsp;n&nbsp;*&nbsp;sizeof&nbsp;(double));</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(yydd&nbsp;==&nbsp;NULL)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="run.c.html#stop">stop</a>("malloc&nbsp;failure");</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;yyqq&nbsp;=&nbsp;(double&nbsp;*)&nbsp;malloc(n&nbsp;*&nbsp;n&nbsp;*&nbsp;sizeof&nbsp;(double));</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(yyqq&nbsp;==&nbsp;NULL)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="run.c.html#stop">stop</a>("malloc&nbsp;failure");</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;initialize&nbsp;D</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(i&nbsp;=&nbsp;0;&nbsp;i&nbsp;&lt;&nbsp;n;&nbsp;i++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;D(i,&nbsp;i)&nbsp;=&nbsp;p1-&gt;u.tensor-&gt;elem[n&nbsp;*&nbsp;i&nbsp;+&nbsp;i]-&gt;u.d;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(j&nbsp;=&nbsp;i&nbsp;+&nbsp;1;&nbsp;j&nbsp;&lt;&nbsp;n;&nbsp;j++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;D(i,&nbsp;j)&nbsp;=&nbsp;p1-&gt;u.tensor-&gt;elem[n&nbsp;*&nbsp;i&nbsp;+&nbsp;j]-&gt;u.d;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;D(j,&nbsp;i)&nbsp;=&nbsp;p1-&gt;u.tensor-&gt;elem[n&nbsp;*&nbsp;i&nbsp;+&nbsp;j]-&gt;u.d;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;initialize&nbsp;Q</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(i&nbsp;=&nbsp;0;&nbsp;i&nbsp;&lt;&nbsp;n;&nbsp;i++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Q(i,&nbsp;i)&nbsp;=&nbsp;1.0;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(j&nbsp;=&nbsp;i&nbsp;+&nbsp;1;&nbsp;j&nbsp;&lt;&nbsp;n;&nbsp;j++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Q(i,&nbsp;j)&nbsp;=&nbsp;0.0;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Q(j,&nbsp;i)&nbsp;=&nbsp;0.0;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;step&nbsp;up&nbsp;to&nbsp;100&nbsp;times</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(i&nbsp;=&nbsp;0;&nbsp;i&nbsp;&lt;&nbsp;100;&nbsp;i++)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(step()&nbsp;==&nbsp;0)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;break;</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(i&nbsp;==&nbsp;100)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="main.c.html#printstr">printstr</a>("\nnote:&nbsp;eigen&nbsp;did&nbsp;not&nbsp;converge\n");</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;p2&nbsp;=&nbsp;D</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(op&nbsp;==&nbsp;EIGEN&nbsp;||&nbsp;op&nbsp;==&nbsp;EIGENVAL)&nbsp;{</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="stack.c.html#push">push</a>(p1);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="tensor.c.html#copy_tensor">copy_tensor</a>();</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p2&nbsp;=&nbsp;<a href="stack.c.html#pop">pop</a>();</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(i&nbsp;=&nbsp;0;&nbsp;i&nbsp;&lt;&nbsp;n;&nbsp;i++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(j&nbsp;=&nbsp;0;&nbsp;j&nbsp;&lt;&nbsp;n;&nbsp;j++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="bignum.c.html#push_double">push_double</a>(D(i,&nbsp;j));</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p2-&gt;u.tensor-&gt;elem[n&nbsp;*&nbsp;i&nbsp;+&nbsp;j]&nbsp;=&nbsp;<a href="stack.c.html#pop">pop</a>();</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;p3&nbsp;=&nbsp;Q</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(op&nbsp;==&nbsp;EIGEN&nbsp;||&nbsp;op&nbsp;==&nbsp;EIGENVEC)&nbsp;{</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="stack.c.html#push">push</a>(p1);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="tensor.c.html#copy_tensor">copy_tensor</a>();</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p3&nbsp;=&nbsp;<a href="stack.c.html#pop">pop</a>();</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(i&nbsp;=&nbsp;0;&nbsp;i&nbsp;&lt;&nbsp;n;&nbsp;i++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(j&nbsp;=&nbsp;0;&nbsp;j&nbsp;&lt;&nbsp;n;&nbsp;j++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="bignum.c.html#push_double">push_double</a>(Q(i,&nbsp;j));</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p3-&gt;u.tensor-&gt;elem[n&nbsp;*&nbsp;i&nbsp;+&nbsp;j]&nbsp;=&nbsp;<a href="stack.c.html#pop">pop</a>();</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;free&nbsp;working&nbsp;vars</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;free(yydd);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;free(yyqq);</tt><br>
<tt>}</tt><br>
<tt></tt><br>
<tt><b><i>//-----------------------------------------------------------------------------</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Example:&nbsp;p&nbsp;=&nbsp;1,&nbsp;q&nbsp;=&nbsp;3</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;c&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;0&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;s&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;0</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;0&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;1&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;0&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;0</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;G&nbsp;=</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;-s&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;0&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;c&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;0</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;0&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;0&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;0&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;1</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;The&nbsp;effect&nbsp;of&nbsp;multiplying&nbsp;G&nbsp;times&nbsp;A&nbsp;is...</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;row&nbsp;1&nbsp;of&nbsp;A&nbsp;&nbsp;&nbsp;&nbsp;=&nbsp;c&nbsp;(row&nbsp;1&nbsp;of&nbsp;A&nbsp;)&nbsp;+&nbsp;s&nbsp;(row&nbsp;3&nbsp;of&nbsp;A&nbsp;)</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n+1&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;row&nbsp;3&nbsp;of&nbsp;A&nbsp;&nbsp;&nbsp;&nbsp;=&nbsp;c&nbsp;(row&nbsp;3&nbsp;of&nbsp;A&nbsp;)&nbsp;-&nbsp;s&nbsp;(row&nbsp;1&nbsp;of&nbsp;A&nbsp;)</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n+1&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;In&nbsp;terms&nbsp;of&nbsp;components&nbsp;the&nbsp;overall&nbsp;effect&nbsp;is...</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;row&nbsp;1&nbsp;=&nbsp;c&nbsp;row&nbsp;1&nbsp;+&nbsp;s&nbsp;row&nbsp;3</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[1,1]&nbsp;=&nbsp;c&nbsp;A[1,1]&nbsp;+&nbsp;s&nbsp;A[3,1]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[1,2]&nbsp;=&nbsp;c&nbsp;A[1,2]&nbsp;+&nbsp;s&nbsp;A[3,2]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[1,3]&nbsp;=&nbsp;c&nbsp;A[1,3]&nbsp;+&nbsp;s&nbsp;A[3,3]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[1,4]&nbsp;=&nbsp;c&nbsp;A[1,4]&nbsp;+&nbsp;s&nbsp;A[3,4]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;row&nbsp;3&nbsp;=&nbsp;c&nbsp;row&nbsp;3&nbsp;-&nbsp;s&nbsp;row&nbsp;1</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[3,1]&nbsp;=&nbsp;c&nbsp;A[3,1]&nbsp;-&nbsp;s&nbsp;A[1,1]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[3,2]&nbsp;=&nbsp;c&nbsp;A[3,2]&nbsp;-&nbsp;s&nbsp;A[1,2]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[3,3]&nbsp;=&nbsp;c&nbsp;A[3,3]&nbsp;-&nbsp;s&nbsp;A[1,3]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[3,4]&nbsp;=&nbsp;c&nbsp;A[3,4]&nbsp;-&nbsp;s&nbsp;A[1,4]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;T</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;The&nbsp;effect&nbsp;of&nbsp;multiplying&nbsp;A&nbsp;times&nbsp;G&nbsp;&nbsp;is...</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;col&nbsp;1&nbsp;of&nbsp;A&nbsp;&nbsp;&nbsp;&nbsp;=&nbsp;c&nbsp;(col&nbsp;1&nbsp;of&nbsp;A&nbsp;)&nbsp;+&nbsp;s&nbsp;(col&nbsp;3&nbsp;of&nbsp;A&nbsp;)</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n+1&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;col&nbsp;3&nbsp;of&nbsp;A&nbsp;&nbsp;&nbsp;&nbsp;=&nbsp;c&nbsp;(col&nbsp;3&nbsp;of&nbsp;A&nbsp;)&nbsp;-&nbsp;s&nbsp;(col&nbsp;1&nbsp;of&nbsp;A&nbsp;)</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n+1&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;n</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;In&nbsp;terms&nbsp;of&nbsp;components&nbsp;the&nbsp;overall&nbsp;effect&nbsp;is...</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;col&nbsp;1&nbsp;=&nbsp;c&nbsp;col&nbsp;1&nbsp;+&nbsp;s&nbsp;col&nbsp;3</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[1,1]&nbsp;=&nbsp;c&nbsp;A[1,1]&nbsp;+&nbsp;s&nbsp;A[1,3]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[2,1]&nbsp;=&nbsp;c&nbsp;A[2,1]&nbsp;+&nbsp;s&nbsp;A[2,3]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[3,1]&nbsp;=&nbsp;c&nbsp;A[3,1]&nbsp;+&nbsp;s&nbsp;A[3,3]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[4,1]&nbsp;=&nbsp;c&nbsp;A[4,1]&nbsp;+&nbsp;s&nbsp;A[4,3]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;col&nbsp;3&nbsp;=&nbsp;c&nbsp;col&nbsp;3&nbsp;-&nbsp;s&nbsp;col&nbsp;1</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[1,3]&nbsp;=&nbsp;c&nbsp;A[1,3]&nbsp;-&nbsp;s&nbsp;A[1,1]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[2,3]&nbsp;=&nbsp;c&nbsp;A[2,3]&nbsp;-&nbsp;s&nbsp;A[2,1]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[3,3]&nbsp;=&nbsp;c&nbsp;A[3,3]&nbsp;-&nbsp;s&nbsp;A[3,1]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[4,3]&nbsp;=&nbsp;c&nbsp;A[4,3]&nbsp;-&nbsp;s&nbsp;A[4,1]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;What&nbsp;we&nbsp;want&nbsp;to&nbsp;do&nbsp;is&nbsp;just&nbsp;compute&nbsp;the&nbsp;upper&nbsp;triangle&nbsp;of&nbsp;A&nbsp;since&nbsp;we</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;know&nbsp;the&nbsp;lower&nbsp;triangle&nbsp;is&nbsp;identical.</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;In&nbsp;other&nbsp;words,&nbsp;we&nbsp;just&nbsp;want&nbsp;to&nbsp;update&nbsp;components&nbsp;A[i,j]&nbsp;where&nbsp;i&nbsp;&lt;&nbsp;j.</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//-----------------------------------------------------------------------------</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Example:&nbsp;p&nbsp;=&nbsp;2,&nbsp;q&nbsp;=&nbsp;5</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;q</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;j=1&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;j=2&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;j=3&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;j=4&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;j=5&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;j=6</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;i=1&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[1,2]&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[1,5]&nbsp;&nbsp;.</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;i=2&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[2,1]&nbsp;&nbsp;A[2,2]&nbsp;&nbsp;A[2,3]&nbsp;&nbsp;A[2,4]&nbsp;&nbsp;A[2,5]&nbsp;&nbsp;A[2,6]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;i=3&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[3,2]&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[3,5]&nbsp;&nbsp;.</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;i=4&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[4,2]&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[4,5]&nbsp;&nbsp;.</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;q&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;i=5&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[5,1]&nbsp;&nbsp;A[5,2]&nbsp;&nbsp;A[5,3]&nbsp;&nbsp;A[5,4]&nbsp;&nbsp;A[5,5]&nbsp;&nbsp;A[5,6]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;i=6&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[6,2]&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[6,5]&nbsp;&nbsp;.</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//-----------------------------------------------------------------------------</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;This&nbsp;is&nbsp;what&nbsp;B&nbsp;=&nbsp;GA&nbsp;does:</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;row&nbsp;2&nbsp;=&nbsp;c&nbsp;row&nbsp;2&nbsp;+&nbsp;s&nbsp;row&nbsp;5</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,1]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,1]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,1]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,2]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,2]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,2]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,3]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,3]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,3]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,4]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,4]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,4]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,5]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,6]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,6]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,6]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;row&nbsp;5&nbsp;=&nbsp;c&nbsp;row&nbsp;5&nbsp;-&nbsp;s&nbsp;row&nbsp;2</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[5,1]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[5,1]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[2,1]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[5,2]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[5,2]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[2,2]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[5,3]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[5,3]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[2,3]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[5,4]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[5,4]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[2,4]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[5,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[5,5]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[2,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[5,6]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[5,6]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[2,6]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;T</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;This&nbsp;is&nbsp;what&nbsp;BG&nbsp;&nbsp;does:</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;col&nbsp;2&nbsp;=&nbsp;c&nbsp;col&nbsp;2&nbsp;+&nbsp;s&nbsp;col&nbsp;5</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[1,2]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[1,2]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[1,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,2]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,2]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[2,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[3,2]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[3,2]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[3,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[4,2]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[4,2]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[4,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[5,2]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[5,2]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[6,2]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[6,2]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[6,5]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;col&nbsp;5&nbsp;=&nbsp;c&nbsp;col&nbsp;5&nbsp;-&nbsp;s&nbsp;col&nbsp;2</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[1,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[1,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[1,2]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[2,2]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[3,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[3,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[3,2]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[4,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[4,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[4,2]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[5,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[5,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[5,2]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[6,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[6,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[6,2]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//-----------------------------------------------------------------------------</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Step&nbsp;1:&nbsp;Just&nbsp;do&nbsp;upper&nbsp;triangle&nbsp;(i&nbsp;&lt;&nbsp;j),&nbsp;B[2,5]&nbsp;=&nbsp;0</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[1,2]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[1,2]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[1,5]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,3]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,3]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,3]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,4]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,4]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,4]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,6]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,6]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,6]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[1,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[1,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[1,2]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[3,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[3,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[3,2]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[4,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[4,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[4,2]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[5,6]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[5,6]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[2,6]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//-----------------------------------------------------------------------------</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Step&nbsp;2:&nbsp;Transpose&nbsp;where&nbsp;i&nbsp;&gt;&nbsp;j&nbsp;since&nbsp;A[i,j]&nbsp;==&nbsp;A[j,i]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[1,2]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[1,2]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[1,5]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,3]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,3]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[3,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,4]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,4]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[4,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[2,6]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,6]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,6]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[1,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[1,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[1,2]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[3,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[3,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[2,3]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[4,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[4,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[2,4]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;B[5,6]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[5,6]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[2,6]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//-----------------------------------------------------------------------------</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Step&nbsp;3:&nbsp;Same&nbsp;as&nbsp;above&nbsp;except&nbsp;reorder</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;k&nbsp;&lt;&nbsp;p&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;(k&nbsp;=&nbsp;1)</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[1,2]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[1,2]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[1,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[1,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[1,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[1,2]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;p&nbsp;&lt;&nbsp;k&nbsp;&lt;&nbsp;q&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;(k&nbsp;=&nbsp;3..4)</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[2,3]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,3]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[3,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[3,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[3,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[2,3]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[2,4]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,4]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[4,5]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[4,5]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[4,5]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[2,4]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;q&nbsp;&lt;&nbsp;k&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;(k&nbsp;=&nbsp;6)</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[2,6]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[2,6]&nbsp;+&nbsp;s&nbsp;*&nbsp;A[5,6]</i></b></tt><br>
<tt><b><i>//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A[5,6]&nbsp;=&nbsp;c&nbsp;*&nbsp;A[5,6]&nbsp;-&nbsp;s&nbsp;*&nbsp;A[2,6]</i></b></tt><br>
<tt><b><i>//</i></b></tt><br>
<tt><b><i>//-----------------------------------------------------------------------------</i></b></tt><br>
<tt></tt><br>
<tt>static&nbsp;int</tt><br>
<tt>step(void)</tt><br>
<tt>{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;int&nbsp;count,&nbsp;i,&nbsp;j;</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;count&nbsp;=&nbsp;0;</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;for&nbsp;each&nbsp;upper&nbsp;triangle&nbsp;"off-diagonal"&nbsp;component&nbsp;do&nbsp;step2</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(i&nbsp;=&nbsp;0;&nbsp;i&nbsp;&lt;&nbsp;n&nbsp;-&nbsp;1;&nbsp;i++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(j&nbsp;=&nbsp;i&nbsp;+&nbsp;1;&nbsp;j&nbsp;&lt;&nbsp;n;&nbsp;j++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(D(i,&nbsp;j)&nbsp;!=&nbsp;0.0)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;step2(i,&nbsp;j);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;count++;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;return&nbsp;count;</tt><br>
<tt>}</tt><br>
<tt></tt><br>
<tt>static&nbsp;void</tt><br>
<tt>step2(int&nbsp;p,&nbsp;int&nbsp;q)</tt><br>
<tt>{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;int&nbsp;k;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;double&nbsp;t,&nbsp;theta;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;double&nbsp;c,&nbsp;cc,&nbsp;s,&nbsp;ss;</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;compute&nbsp;c&nbsp;and&nbsp;s</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;from&nbsp;Numerical&nbsp;Recipes&nbsp;(except&nbsp;they&nbsp;have&nbsp;a_qq&nbsp;-&nbsp;a_pp)</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;theta&nbsp;=&nbsp;0.5&nbsp;*&nbsp;(D(p,&nbsp;p)&nbsp;-&nbsp;D(q,&nbsp;q))&nbsp;/&nbsp;D(p,&nbsp;q);</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;t&nbsp;=&nbsp;1.0&nbsp;/&nbsp;(fabs(theta)&nbsp;+&nbsp;sqrt(theta&nbsp;*&nbsp;theta&nbsp;+&nbsp;1.0));</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if&nbsp;(theta&nbsp;&lt;&nbsp;0.0)</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;t&nbsp;=&nbsp;-t;</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;c&nbsp;=&nbsp;1.0&nbsp;/&nbsp;sqrt(t&nbsp;*&nbsp;t&nbsp;+&nbsp;1.0);</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;s&nbsp;=&nbsp;t&nbsp;*&nbsp;c;</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;D&nbsp;=&nbsp;GD</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;which&nbsp;means&nbsp;"add&nbsp;rows"</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(k&nbsp;=&nbsp;0;&nbsp;k&nbsp;&lt;&nbsp;n;&nbsp;k++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;cc&nbsp;=&nbsp;D(p,&nbsp;k);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;ss&nbsp;=&nbsp;D(q,&nbsp;k);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;D(p,&nbsp;k)&nbsp;=&nbsp;c&nbsp;*&nbsp;cc&nbsp;+&nbsp;s&nbsp;*&nbsp;ss;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;D(q,&nbsp;k)&nbsp;=&nbsp;c&nbsp;*&nbsp;ss&nbsp;-&nbsp;s&nbsp;*&nbsp;cc;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;D&nbsp;=&nbsp;D&nbsp;transpose(G)</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;which&nbsp;means&nbsp;"add&nbsp;columns"</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(k&nbsp;=&nbsp;0;&nbsp;k&nbsp;&lt;&nbsp;n;&nbsp;k++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;cc&nbsp;=&nbsp;D(k,&nbsp;p);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;ss&nbsp;=&nbsp;D(k,&nbsp;q);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;D(k,&nbsp;p)&nbsp;=&nbsp;c&nbsp;*&nbsp;cc&nbsp;+&nbsp;s&nbsp;*&nbsp;ss;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;D(k,&nbsp;q)&nbsp;=&nbsp;c&nbsp;*&nbsp;ss&nbsp;-&nbsp;s&nbsp;*&nbsp;cc;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;Q&nbsp;=&nbsp;GQ</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<b><i>//&nbsp;which&nbsp;means&nbsp;"add&nbsp;rows"</i></b></tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for&nbsp;(k&nbsp;=&nbsp;0;&nbsp;k&nbsp;&lt;&nbsp;n;&nbsp;k++)&nbsp;{</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;cc&nbsp;=&nbsp;Q(p,&nbsp;k);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;ss&nbsp;=&nbsp;Q(q,&nbsp;k);</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Q(p,&nbsp;k)&nbsp;=&nbsp;c&nbsp;*&nbsp;cc&nbsp;+&nbsp;s&nbsp;*&nbsp;ss;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Q(q,&nbsp;k)&nbsp;=&nbsp;c&nbsp;*&nbsp;ss&nbsp;-&nbsp;s&nbsp;*&nbsp;cc;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;}</tt><br>
<tt></tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;D(p,&nbsp;q)&nbsp;=&nbsp;0.0;</tt><br>
<tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;D(q,&nbsp;p)&nbsp;=&nbsp;0.0;</tt><br>
<tt>}</tt><br>
<tt></tt><br>
</body></html>
